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Thermodynamical properties of the helix-coil transition were successfully described in the past 
by the model of Lifson, Poland and Sheraga. Here we compute the corresponding structure factor 
and show that it possesses a universal scaling behavior near the transition point, even when the 
transition is of first order. Moreover, we introduce a dynamical version of this model, that we 
solve numerically. A Langevin equation is also proposed to describe the dynamics of the density 
of hydrogen bonds. Analytical solution of this equation shows dynamical scaling near the critical 
temperature and predicts a gelation phenomenon above the critical temperature. In the case when 
comparison of the two dynamical approaches is possible, the predictions of our phenomenological 
theory agree with the results of the Monte Carlo simulations. 

PACS: 87.10.Gg, 64.60.Fr, 64.60.Ht 



I. INTRODUCTION 

Ever since the discovery of the helicoidal structure of 
DNA under the usual conditions, the phenomenon of its 
transformation to a coil under heating has been studied 
and described as an equilibrium phase transition. Basic 
works by Lifson and subsequently by Poland and Sheraga 
led to a model (LPS) considered today standard [J 
The model is one-dimensional and shows a phase tran- 
sition only because the interactions are long-ranged. A 
more realistic description of the two strands moving and 
interacting in three dimensional space has also been pro- 
posed and studied numerically [2 [H, ] • A model of 
LPS-type was shown to be able to consistently fit all rel- 
evant simulation results Q . More recently, some authors 
have proposed a new mechanism which yields an abrupt 
first-order phase transition [1, at variance with the 
previous theoretical considerations in favor of a continu- 
ous one. On the experimental side the situation is still 
far from conclusive. Previous measurements indicated 
the presence of a sharp jump in the fraction of bound 
base pairs [13, EH ■ More recently, the results reported in 
[l2| have been interpreted as an indication of a weakly 
continuous phase transition. 

It seems, however, that little is known about the dy- 
namics of this denaturation transition. We can mention 
the contributions based on a mechanical approach [T3l.[l4j 
and on a Langevin description of the two DNA-strands 
as polymers in continuous space [TH - 

In this paper we investigate the original LPS model 
while considering the exponent a as a free parameter. 
For what concerns the equilibrium properties we compute 



the structure factor and find that near the critical point it 
exhibits universal scaling behavior, which depends only 
on a. Moreover, we also look at the transition from a 
dynamical point of view. We do it in two ways: first we 
describe this dynamics by a master equation such that 
the equilibrium state is chosen to be the one of the LPS 
model. Secondly, we introduce a Langevin type equation 
for the order parameter which is the density of hydrogen 
bonds. The analysis of the solution of this equation re- 
veals a critical slowing-down, the relaxation time diverg- 
ing at the critical point with a power law still depending 
on a. In the absence of noise, we find that above the crit- 
ical temperature a gelation phenomenon occurs, namely 
the density of H-bonds approaches zero in a finite time. 
With the addition of noise we obtain also the analytic 
expression for the power-spectrum of the density of H- 
bonds. Finally, we compare the two descriptions and find 
that the dynamics of the order parameter derived from 
the master equation and from the Langevin equations 
essentially agree. 

The following three sections are devoted to a study of 
equilibrium properties. In Section II we summarize the 
basics of the LPS model. Formulas for the free energy 
and the pressure are given here, as well as scaling formu- 
las for the density in the vicinity of the phase transition. 
Section III contains the computation of the one- and two- 
point correlation functions for finite systems and in the 
thermodynamic limit. Our results for the structure fac- 
tor are presented in Section IV. Results on the dynamics 
are collected in Sections V and VI. In Section V we de- 
scribe the microscopic dynamics realized numerically in 
this work. The equation for the order parameter and its 



2 



analysis both with and without noise is given in Section 
VI, where we also compare the results of the two dynam- 
ical approaches. The paper ends with a Summary. 



II. THE LPS MODEL 

Let and 1 denote respectively a broken or existing 
H-bond at a given site. Then, a configuration of the 
double-strand of DNA of length TV is represented by a 
sequence n = {n x } x=1 where n x = or 1. In n there are 
alternating maximal connected subsequences (intervals) 
of 1 and 0, respectively. A weight a(j) > or b(j) > is 
associated to an interval of length j of 1 or 0, respectively. 
If the lengths of the intervals of n are k and rrij then the 
weight of n is 



and 



-V(n) 



Hamlin? 



(2.1) 



As usual, thermal equilibrium can be discussed in the 
canonical or in the grand-canonical ensemble. In the 
first case the description is restricted to configurations 
with the total number of H-bonds, n — ^ x n x , fixed. In 
the second case all configurations compatible with the 
boundary condition are considered but an activity 
is introduced. The chemical potential p is the energy 
needed to break a H-bond. The canonical partition func- 
tion is 



IN.' 



-V(n) 



(2.2) 



and the grand-canonical partition function is 

Q N (m =J2QN, n e n ^ = J2e- u ^ (2.3) 

n n 

where C/(n) = V(n) — /3p^2 x n x . Alternately, e~ u ^ is 
obtained from e~ v ( n ^ by replacing each factor a(l) by 
a(l)e 1 ^^. The free energy / and the pressure p are re- 
spectively given by 



and 



/?/ = - lim — IiiQn n 

N -> oo N 
n/N -> p 



f3p= lim — lnQ N (Pii). 

N—>oc IV 



(2.4) 



(2.5) 



The temperature does not appear explicitly in the canon- 
ical partition function, therefore /?/ is independent of f3 
and Pp depends on it through e = Pp. The free energy 
and the pressure are Legendre transforms of each other, 

pp(e) = sup{ep - Pf(p)} = e Pe - 0f(p e ) (2.6) 
p 



Pf(p) = sup{ep - (3p(e)} = e p p - (3p(e f 



(2.7) 



Here p is the density of 1 (H-bonds), and p e and e p are 
the respective values at which the suprema are attained. 

The partition functions depend on the boundary con- 
ditions imposed on the system. Typical choices are the 
free and the periodic conditions, or fixed or 1 on the 
boundary sites. For example, with the choice ... 1 or 
1...0, 



Q N ,n = 22M n ) B k(N -n) 



k>l 



where 



and 



h>l l k >li=l 



(2. 



(2.9) 



B k (n)= ■■■ J2 ]jHmi)S n ^ mi - (2.10) 

mi > 1 nih > 1 i—1 

Ak{n) = Bfc(n) = if n < or k > n, therefore Qn,u — 
for n > N. One can show that in the thermodynamic 
limit the free energy and the pressure do not depend on 
the boundary condition. Below we give these quantities 
in the case of the simplest model (the LPS model in a 
strict sense) obtained by choosing 

a(l) = a b(l) = bu l r a (2.11) 

where u, a > 1 and a, b > . With the notations c = ab, 



and 



one obtains 



Inu-ln[l + cCa(0)], 



(3p = In ti + 0(e — e c )x(e) 



(2.12) 



(2.13) 



(2.14) 



where O is the Heaviside function and, for e > e c , a;(e) is 
the unique solution of the equation 



1 + c ( a (x) . 



(2.15) 



Thus, x(e c ) = and x(e) is monotone increasing with 
e. For e < e c , (3p takes on its minimum value, In it. 
Comparison with ([2^4]) . ([2~5| and (|2.11|) shows that the 
corresponding density and free energy are p e = and 
f(p e ) = — hiu. Legendre transformation yields f(p) for 
all p. With p c = p ec , 



(3f( P ) = 



-lnu+e c p, p < p c 

-]nu-y + p(y + \n 1+c u Ca{v) ) , P>Pc 

(2.16) 
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where y = y(p) is the solution of the equation 
, , cCa-i(y) _ „_ a 



l + cCa(y) 



(2.17) 



In fact, y(p) = x(e p ); y(p) is defined for p c < p < 1 and 
is monotone increasing, y(p c ) — and y(l) = oo. There- 
fore, f{p) is monotone increasing from /(0) = — f3 lnii 
to /(l) = 0; as the Legendre transform of the convex 
p(e), f(p) is also convex. Substituting y = into (|2.17|) 
yields 



l + cCc(0) 



l + c[C o -l(0) + Ca(0)]' 



(2.18) 



This shows that p c = if a < 2 and p c > if a > 2. The 
first case corresponds to a continuous phase transition 
with the order parameter p e vanishing continuously as e 
approaches e c from above. In the second case the transi- 
tion is of first order: p e decreases continuously to p c as e 
decreases to e c and jumps to zero below e c . Computation 
of p' e — dp € / de from Eq. (|2.17p shows that p' ec+Q = oo if 
a < 3 (while it is finite for a > 3). The type of divergence 
can be extracted from the asymptotic forms 



P< 



-(e-e c ) 
ln(e- 



2-a 
2 



if 1 < a < 2 

if 2 < a < 3 
if a = 2 



(2.19) 



valid when e decreases to e c . Equation (|2.19[) is obtained 
from the asymptotic form of ( a (x) for small x and from 
x(e) w e - e c for e - e c < 1. It will be used in Section 
Vll Finally, we note that for e < e c (T > T c ), when 
/9 £ = 0, the decay of the density with increasing N can 
be computed, 



Pn 



= (l-T c /T))jV 



(2.20) 



This formula is valid for N > (1 - T c /T)" 



III. CORRELATION FUNCTIONS 

We will compute the one- and two-point correlation 
functions of H-bonds in the grand-canonical ensemble for 
the LPS model. Let us denote the partition function of 
the box [x, y] by 



Q a ^x lV ] = Q aP (v-x + l) 



(3.21) 



where we have imposed the boundary conditions a at i 
and at y (Note that Q al3 = Q' 3 ", Q 00 (l) = bu and 
Q 11 (l) = ae e ). Without prejudice of generality we com- 
pute the correlation functions with the boundary condi- 
tions a = 1 and (3—1. We begin with the density Pn{t)- 
In a configuration the site r is in a cluster of l's begin- 
ning at a site x + 1, ending at a site y — 1 , so that at sites 



x and y we have 0's. It follows that 

JV-l r-l 

Q 11 {N)p N {r)=a £ £ Q 10 [l, xje^—^g 01 ^ tf] 

y—r+l x— 2 



By introducing 



,(iV) = ^Q 01 (x)e- 



x=2 



we can rewrite (I3.22[) as follows: 

ae eN [g(r - 1) + %(iV -r) + 1] 



PJv(r-) 



(3.22) 
(3.23) 

(3.24) 



For the two-point correlation function p^{r,r') we ob- 
serve that in a configuration either r and r' are in the 
same cluster of l's or they belong to distinct clusters of 
l's. By assuming r' > r + 1 we obtain 



Pn(t, r') 



ae eN [g{r - 1) + 1] [g(N - r') + 1] 
WW) > 

r' — r 

l + a^2 Q 00 (l)e~ d (r' -r-l) 
i=i 



(3.25) 



These expressions can be simplified by considering that 
the following relations hold for N > 2: 



(JV) = [u 2 e- 2e ZAr +1 - 2ue- e ZAT + Z N -t] 



(N) = au N Z N 



where 



Z 



1 



N 



dw 



1 



2iri J c w N+1 ue~ e - R(w) 
C is some circle around the origin and 

oo 

R(w) = w + c ^2r a w l+1 
1=1 

Upon these results we have 

g(N) + l = u N+1 e-^ N+ ^Z_ 
and finally 

: Z r —\ZN—r 



N 



p N (r) = ue 



... 9 —2(:Z r —\Z]\[ — r 'Z r '— r 

pN{r,r) = u e 

"N-l 



(3.26) 
(3.27) 

(3.28) 

(3.29) 

(3.30) 
(3.31) 



where r' > r. 

In the low-temperature phase (e > e c ) we can perform 
the thermodynamic limit (N — * oo) by fixing \r — r'\ and 
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letting the distance of r to the boundaries also tend to 
infinity. Making use of the asymptotic (N 3> 1) relation 



with k > 0, 



and 



Z N ~ \ N q(l + e- KN ) 



\- L = -w* 
u 



, dfl(«Q . 
dw; 



(3.32) 



(3.33) 



(3.34) 



where u>* is the positive solution of the equation ue e = 
R(w*), we obtain 

p(r) = p = ue~ e q = d(3p/de, (3.35) 

the expected density of H-bonds, and 

p(r,r') = qu 2 e~ 2t X^ r '- r ^Z lr ,_ rl . (3.36) 

We can see that the truncated correlation function 
p T (r, r') = p(r, r') — p(r)p(r') decays exponentially 



p (r,r') ~ p e 



2 -«j| r -r| 



(3.37) 



for \r' — 



oo. We note that translation invariance 



is lost in the high-temperature phase, because there re- 
mains an explicit dependence on r and r 1 in p(r, r'), even 
if r and r' are far from the boundaries. This can be con- 
sidered as a weakness of the LPS model. Nonetheless, 
we shall see that we can still attribute a clear meaning 
to the structure factor in this phase. 



IV. STRUCTURE FACTORS 

The structure factor is more accessible to experiments 
than the corresponding correlation function. It is defined 
as follows 



(4.38) 



We aim at computing the limit S(k) of this quantity as 
N — > oo. The two cases T <T C and T > T c are treated 
separately. 

Below T c — The correlation functions become transla- 
tional invariant for N — > oo, i.e. p]^(r,r') p T (0,r' — r). 
Moreover p T (0, r) decays exponentially and we obtain 



Making use of (pT2Tj) . (|3~2"8) and (|3~37|) we have 



S(k) 



2Re 



l + c( a (6)-e* k [l + c<; a (S-ik)} 



(4.39) 



- 1. 

(4.40) 



In this expression we have reparametrized w* as 

w*=e~ & . (4.41) 

Accordingly, T — T c corresponds to 6 = 0. The ex- 
pression (|4.40[) depends on the details of the LPS model. 
However, we may expect that in the critical region, char- 
acterized by |fc| << 1 and | T y T ° | << 1 universal features 
emerge and the results depend only on the parameter a . 
If this is the case, we can be confident of the predictions 
of the model. Let us focus on the analysis of the critical 
region. The parameter S allows us to define a correlation 
length £ through the relation 



Near T c , we find 



£o(a)(T c /T-l) 



if 1 < a < 2 



£i(a)(T c /T- I)" 1 if2<a<3. 



(4.42) 



(4.43) 



In order to obtain explicit expressions we need to study 
the function ( a (5), cf. Eq. (|2.12p . for small values of S. 
For this purpose we can use the integral representation 



US) = 



i 



t 



a-l 



■dt 



r(a) J e t+A - 1 
valid for Re S > 0. By this relation we can write 



(4.44) 



Ca{8) - C*(8 - *k) (4-45) 
= (e- lk - 1) [Ca-i(S) + (l-e~ lk )B a (k,6)] 



where 
B a (k,S) 



T(a) 



■dt. 



For 1 < a < 3 and S 
to 



(4.46) 

the above expression reduces 



B a {k,8) 



?Q-3 



with 



b a (k) 



T(2-a) 



(i-iky 



b a (k) 



- 1 



(4.47) 



In particular, 



MO) 



A- 2 



P(3 



+ i 



.r(2-a) 



(4.48) 



(4.49) 



The compressibility x = S(0) can now be computed. For 
1 < a < 2 the phase transition is continuous and we have 



P 

X 



PoC- 2 
Po(2 - a)£ 



2a-3 



(4.50) 
(4-51) 



For 2 < a < 3 the phase transition is of first order and, 
as p — > p c , one has 



X ~ Xo£ 



3-a 



(4.52) 
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We can also find the scaling limit for S(k). It is denned 
by taking the limits k — ► and £ — » oo with the product 
/c£ fixed at some finite value. If 1 < a < 2 we find 



5(fc) 2(a-l) 



Re 



1 



(1 - iki) 01 - 1 - 1 



5(0) 2- a 

For small fc£ this gives a Lorentzian behavior 

5(fc) _ 4 
5(0) ~ (feO 2 + 4 

while in the critical region one has 

S(k) _ 2{a- l)cosf (a- 1) 
S(0) ~ (2-a)(A;0 Q - 1 ' 

On the other hand, for 2 < a < 3 we have 
S(k) _ 2 Re [1 - (1 - ^C)"" 1 ] 



5(0) 



(a-l)(a-2)(A;0 s 



(4.53) 



(4.54) 



(4.55) 



(4.56) 



For |fc£| << 1 one has again a Lorentzian, but in the 
critical region (|4.56[) takes the form 



S(k) _ 2 cosf (a+ 1) 
5(07 ~ (a-l)(a-2)(fe0 3_Q ' 



(4.57) 



Above T c - - In this phase the properties of the 
LPS model are more dependent on N than in the low- 
temperature one. Basically, one has to provide a suit- 
able description of the partition function Zn, defined in 
(|3.27[) . Let us introduce the parameter v = ue~ e . It can 
easily be shown that if v > v c then 



Z7T 



+ 7V 



-iN8 



where 



0-1(0) =u _e i9 [l + cC a H6 



(4.58) 



(4.59) 



In this case the structure factor can be expressed in the 
form 

NS N {k) = 2v 2 ReC N (k)-v 2 \D N {k)\ 2 ~vD N {0) (4.60) 
where 



C N (k) = Z N \ 
and 

D N (k) = Z N \ 



+ 7T 



d0 
2T 



-i(N-l)6 2 



g 2 (8)g(6 + k) (4.61) 



+ 7T J/J 

^ e -^-^ 9 (0) ff (0 + fc). (4.62) 



In order to find the asymptotic behavior (N — > oo) of 
these expressions we use the following formulae. Let F(x) 
be a function whose only singularity in the interval [—a, a] 



is at the origin and the singularity is integrable. If, close 
to the origin, 



then 

/ F(x)e- iNx dx ~ ^^(s)r(l + r J )(iVe^ s )-( 1+r ^ , 

(4.64) 

where r*j is a sequence of real numbers such that — 1 < 
r\ < r 2 < • • • . Using Eqs. (|4.44j) - (|4.48j) . this formula 
yields 

Z N ~ c/{S 2 N a ) (4.65) 
where 5 = v — v c . Moreover, if k ^ then 

C N (k) ~ fl (fc) + ff 2 (-fc)e ife ( JV - 1 ' (4.66) 

and 

D N (k) ~ g ( k ) + gi-ky^-V. (4.67) 
This shows that the compressibility 

X = 5(0) - 2vv c /{NS 3 ). (4.68) 
Accordingly, in the limit N — > oo we have 

§^ = ^| ff (fc)| 2 [(2/5)Re 5 - 1 (fc) - 1] - A. (4.69) 
5(0) w c w c 

Since 5 — 8q(T/T c — 1), also in this case we can obtain the 
scaling limit. With the definition (|4.43|) of the correlation 
length £, for 1 < a < 2 the normalized structure factor 
reads 



S(k) 



l + 2(fc£)"~ 1 sinfa 



5(0) l + (^) 2 ( Q - 1 )+2(fc^)«- 1 sinfa 

while, for 2 < a < 3, 

5(fc) _ 1 
5(0) " T+W 



(4.70) 



(4.71) 



We stress that the main result of these computations is 
the existence of a universal scaling limit for the normal- 
ized structure factor 5(fc)/5(0), which depends only on 
the parameter a. This suggests also that the LPS model 
belongs to a new universality class, in the language of 
critical phenomena. Accordingly, its predictions about 
the process of DNA denaturation can be taken with some 
confidence. Moreover, we expect that this result may 
help to clarify experimentally the nature of the phase 
transition, cither continuous or first-order. 
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V. MICROSCOPIC DYNAMICS 

We will describe the evolution of the configuration n 
by a discrete-time local dynamics. Although the method 
is well-known, we present it in some detail. 

In a time step, n can change only at a single site x. 
The new configuration is T x n and, hence, (T x ii) y = n y 
if y 7^ x, and (T x n) x = 1 — n x . Note that T x equals 
the identity. The transition probability from n to m is 
Pn.m- More precisely, it is the conditional probability 
that, if the configuration is n in time t, it will be m in 
time t + 1. It then follows that p n .m = unless m = n 
or m = T x ii for some x. The irreducibility of the matrix 
Pn.m and, hence, the ergodicity of the dynamics and the 
uniqueness of the equilibrium state is guaranteed if p n ,T w n 
is indeed positive for all x. (Note, however, that both 
n and T x n have to satisfy the boundary conditions. If 
these are given by fixing ri\ and tin then 2 < x < N — 1.) 
In a numerical experiment the transition probability is 
decomposed as 



Pn,r in = p x W(T x n\n, x) 



(5.72) 



where p x is the probability to choose x for an attempt 
of change, and ^(T^nln, x) is the conditional probabil- 
ity that once n and x have been selected, the change is 
executed. Clearly, J2 X P X = 1 m ust hold; for periodic or 
free boundary conditions the typical choice for p x is 1/N. 
The probability to stay in n if the change is attempted 
in x is 



W(n\n, x) = 1 - W(T x ii\n, x) 



and, thus, 



N 



N 



(5.73) 



(5.74) 



x=l 



This is nonzero unless PL(n|n,:E) = for all x. 

Let pt(n), t = 0, 1,2, . . ., be the probability of n at 
time t. The master equation of the evolution is 



A' 



Pt+i(n) =Pt(n)p„.„ + y^p t (r x n)p Txn;n 



(5.75) 



x=l 



The local transition probability W is to be chosen so 
that the grand-canonical distribution n be the unique 
equilibrium solution of (|5.75p . Imposing the condition of 
detailed balance 

7r(n)W(T x ii\n, x) = n(T x n)W(n\T x n, x) , (5.76) 

we still have an infinity of choices: any W satisfying 

W(T x n\n,x) 



W(n\T x ii,x) 



z(n,x) (5.77) 



and < W^T^nln, x) < 1 for all x and all n of length N 
with the right boundary conditions will do. 



If 7 is an interval (a sequence of consecutive integers 
between 1 and N) then by definition the distance of x 
to 7 is d(x,j) — min{|x — y\ : y € 7}. Because U(n) is 
additive in the contributions of the intervals of n, 



U(n) = - ^2[bxa(li)+lie] - ^ln6(m 3 ) 



(5.78) 



z(n, x) depends only on the intervals of n and T x ii whose 
distance to x is or 1 . If in n there is a single interval 
with this property, there will be three in T x ii, and vice 
versa; if there are two in n, there will also be two in 
T x n. We obtain altogether three different functions of 
the lengths of intervals, and their reciprocals. These are 
listed in Table 1. 



(1)^1(1)^ -> (l)'iO(l)'» 


Ri 


(0) h l(l) 12 -> (0) il 0(l)' 2 


R2 


(l) ;i l(0) !2 -> (l) !l 0(0) i2 


R 3 


(0)' 1 1(0)' 2 -> (0) !l 0(0)' 2 


Ra 


(l) ! i0(l) i2 - (1)^1(1) !2 




(0)' 1 0(1)' 2 -> (0) !l l(l)' 2 


R 2 


(i)^o(o) i2 -» (i) l n(iy* 




(0) !l 0(0)' 2 -» (0) !l l(0) !2 


RI 1 



TABLE I: Transition rates z(n,x). The notations {Xf 1 and 
(A)' 2 indicate a left-cluster of length h and a right-cluster 
of length I2 made of A-symbols, respectively. R\ — abv, 
R 3 = "(77tT) a , ^ = ^ (j^T where 



R2 = v ( 
v = ue~ e 



h y 



Typical choices of W are 

W(T x n\n,x) =mm{l.z(n,x)} or ^"'^ N . (5.79) 

1 + z(n, x) 



They satisfy (|5.77|) because z(T x n,x) — l/z(n, x). We 
shall work with the first one. 



VI. THE ORDER PARAMETER EQUATION 

The microscopic dynamics of denaturation as de- 
scribed by the master equation is quite complex and un- 
tractable analytically. Most of the time, however, one 
is interested in the evolution of the order parameter, 
namely the density p of 1 (H-bonds). This quantity, ob- 
tained as an average over the sample, should vary much 
more slowly than the other degrees of freedom to which 
it is coupled. This coupling has two effects, (i) It creates 
an effective thermodynamic force exerted on the order 
parameter and depending nonlinear ly on it. This force 
should vanish if the density, and via it the pressure, reach 
their equilibrium value, corresponding to the prescribed 
value of e. We choose it therefore in the form 
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where f(p) is the equilibrium free energy for each value 
of p, cf. Eqs. (|2.16[) . So the deterministic part of the dy- 
namics realizes the search of the supremum in Eq. (|2.6|) . 
One can recognize here also the mean-field approach to 
the problem of the critical slowing down in phase transi- 
tions, (ii) Due to the very large number of the coupled 
degrees of freedom, the coupling creates a random force 
of zero average, responsible for the fluctuations of the 
density. We choose the random force as a white noise 
rj{t) with a correlation 



(r,(t)r,(t')) = 2 7 (3- 1 8(t-t?). 



(6.80) 



Thus, the time evolution of p(t) is governed by the 
Langevin equation 



p = -7[(Af)'(p) - IM + \T$ v(t) 



(6.81) 



where 7 > 0. We expect this equation to describe cor- 
rectly the behavior of the order parameter in the vicinity 
of the transition point. The comparison of the analytical 
predictions based on Eq. (|6.8ip with the numerical solu- 
tion of the master equation can be a first check of this 
hypothesis. From now on we consider only the LPS case 
(HOT]) . From Eq. (|2~TrJ)) we deduce 



W)'(p) = e„, 



therefore (|6.81[) reads 



P(t) = -7k,« - e] + V(t) 



(6.82) 



(6.83) 



As we have seen in the former section, the microscopic 
dynamics depends only on a, c = ab and v. It is im- 
portant to realize that e p — e also depends only on these 
parameters. Indeed, for given p we solve Eq. (|2.17|) for 
y, replace x by y(p) in Eq. (|2. 15|) and extract e p . Thus, 



- e = y(p) - ln{l + cC a [y(p)}} + hi'- 



(6.84) 



On the other hand, we expect to obtain the value of 7 
from the microscopic dynamics. 

Note that e p = e c if p < p c and tends to 00 as p 
approaches 1 . Inverting p e shows that e p starts with zero 
derivative if a < 3, while e Pc+ o > if a > 3. We will 
confine ourselves to the case a < 3. 



A. Evolution without noise 

CASE 1: 2 < a < 3. In this case p c > 0. The evo- 
lution is naturally different above and below the critical 
temperature. 

(i) T > T c (e < e c ). Let p = p{0). We distinguish 
between two subcases. 

1. p < p c . Then p = -7e c (l - T c /T) and 

pit) = po - 7£ C (1 - T c /T)t for 0<t<t f (6.85) 
where tf = po/^e c (l — T c /T) and pitf) = 0. 



2. po > p c - Then for < t < t c , p evolves according to 
p = — 7(e p — e), where t c is defined by p(t c ) = p c . For 
t > t c the force is constant as in 1., hence 

p(t) = p c -je c (l-T c /T)(t-t c ) for t c <t<t f (6.86) 

where tf = t c + p c /je c (l — T c /T) and p(tf) = 0. 

What we see here is the phenomenon of gelation: The 
system reaches equilibrium in a finite time [l6| . 

(ii) T < T c (e > e c ). Then ]im t ->oo p(t) — Pe > Pc- 
For large t, by linearizing e p about e we find e p — e w 
[p — p e )/p' e [recall that p' e = dyO e /de] and 



pit) - Pe - Ae- h/p '- ]t . 



(6.87) 



When T is close to T c , from Eq. (|2.19p we deduce p' e 
ia - 2)e c (T c /T - 1) Q " 3 and 



pit) - p e ~ exp(-t/r) 



(6.88) 



with 



T = p'Jj~J- 1 (a-2)iT c /T-l) 



a-3 



CASE 2: 1< a < 2, p c = 0. 

(i) T > T c . Now starting from p = p(0) > 0, for all 
times e p > e c > e, therefore p < — 7e c (l — T c /T) and 
pit) attains zero in a finite time tf < po/7 £ c(l — T c /T). 
Again, we find gelation. 

(ii) If T < T c , for large times we find the same result as 
in (|5^T|) and (jrTgg)) with 



r = pi(e)h~ 1 - 1 i2- a )iTc/T-l) 



l-a 



CASE 3: a — 2. For T > T c we obtain gelation, for 
T < T c an exponential relaxation to poo = p t with decay 
time 

T ~ [l{T c /T — 1) ln 2 (T c /T — 1)] — 1 . 

As in the case of the static structure factor, we can 
obtain a scaling limit for the time-dependent order pa- 
rameter p{t). It can be determined by taking p — > 0, 
t — » 00, T — > T c . We restrict this analysis to the case 
a > 3/2. 

(i) T > T c . The gelation time i/ diverges as 



tf - (T/T c - iy (1 ~ u) where f 
(ii) T <T C . Then 

^=zte 1 _ T/Tc) i 

where z(t) is the unique solution of 

z = -{z 1/v - 1) 



2- a 
a -I' 



(6.89) 



(6.90) 
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tending to 1 as t tends to infinity. Therefore, if i(l 
TjT c ) l - v > 1 then 



pit) ~ p t 



1 + cxp 



1±_ 



(l-T/T c )t 



(6.91) 



On the other hand, in the critical region 

l«t« \T/T C - 

one has 

-u/(l-u) 



p(t) 



1 - v 



(6.92) 



B. The effect of noise 

In order to compute the time-dependent correlation 
functions of the density p, we have added a white noise 
term rj(t) to the deterministic order parameter equation 
(|6.8ip . On the other hand, this choice may lead to a 
consistency problem. In fact, it is well known that in 
this case the stationary distribution V(p) of the density 
is of the form 



P(p)~exp }--\J3f(p)- £ p] 



(6.93) 



The average of p taken with V(p) is, in general, different 
from the equilibrium value p c , given by the equation 



(Pf)'(pe) = e, 



(6.94) 



cf. Eq. (|2.6|) . This is a direct consequence of having 
interpreted f(p) in (|6.8ip as the equilibrium free energy. 
Therefore, we will only investigate the effect of noise by 
considering the linearized version of Eq. (|6.81[) close to 
equilibrium, i.e. 

P = -l(Pf)"(Pe)(p - Pe) + V%(*). (6.95) 
Below T c , f"(p) is related to the compressibility, 
d 2 f 



dp 2 



k B T/ X 



(6.96) 



which is a well defined quantity in the limit N — > oo. This 
is not the case above T c , where the deterministic equation 
predicts the occurrence of gelation. In order to circum- 
vent this difficulty, we assume that (|6.95|) and (|6.96p are 
still valid above T c , while we maintain the explicit de- 
pendence on N of both p e and x- The time-dependent 
correlation function C p (t) is defined as 

C p (t) = (p(t + t')p(t')) (p(t + t')){p{t')) (6.97) 

where the bracket denotes the time average. The power 
spectrum of p is given by 



1 r + °° 
W = ^ J e wt C p (t)dt 



(6.98) 



Making use of (|6.95p and (I6.96P one obtains 

C P W=xe- |t|/r (6.99) 

and 

XT 1 



7r 1 + [tory 

where r = xll- F° r T /" T c we obtain 

_ f 7~ 1 <f "~ 3 if 1 < a < 2 
T ~ I l~H 3 ~ a if 2 < a < 3 



(6.100) 



(6.101) 



with £ defined in Section IIV1 In the language of critical 
phenomena Eq. (|6.10ip can be interpreted as a scaling 
law valid for lot not too large. The dynamical critical 
exponent z is then given by 



= I 2a - 3 if Ka<2 
Z ~ I 3 - a if 2<a<3. 

Near to but above T c , one has instead 
since x = 2vv c /NS 3 , see f)4.68[) . 



(6.102) 



(6.103) 



C. Comparison with the numerical experiment 

We have performed numerical simulations according to 
the microscopic MonteCarlo dynamics described in Sec- 
tion [Vl in order to check its agreement with the order- 
parameter equation. We report here the results obtained 
for the two cases a — 2.5 (first-order phase transition) 
and a = 1.5 (second-order phase transition). The mi- 
croscopic dynamics is found to depend on the control 
parameter v = ue~ e . Since e = f3p, v can be used at the 
place of the temperature T ~ /3 _1 for describing numer- 
ical results. Without prejudice of generality one can fix 
ab = 1. According to Eq. (|2. 13(1 the critical value v c is 
given by the relation 



v c = 1 + C«(0) 



(6.104) 



Most of the results discussed in this section have been ob- 
tained for lattice size N — 5 x 10 3 . The time evolution of 
p(t) has been extended over time spans ranging between 
lOOr and 600r. In this subsection r = N denotes the 
natural time unit of lattice updates (not to be confused 
with the r used in Eq. (|6.99p ) . Fluctuations have been 
smoothened by averaging over a large number of initial 
conditions (typically 10 4 ). The initial conditions have 
been sampled among high density equilibrium states ob- 
tained for v = 0.1 (p(0) s» 0.9). For the special situation 
concerning the case a = 2.5 and p(0) < p c the initial con- 
ditions have been sampled by attributing to each site the 
value "1" with probability p(0) and "0" with probability 
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(1 — p(0)) ■ Lacking an equilibrium state of finite den- 
sity due to the first-order nature of the phase transition, 
this is a natural choice for sampling states of fixed ini- 
tial density p(0). We have also verified that other recipes 
can modify the duration of the transient evolution, but 
we have not observed any difference in the long time be- 
havior. Finally, we have also verified that the results do 
not depend on the choice of boundary conditions. In par- 
ticular, the results reported in this subsection have been 
obtained for fixed boundary conditions, where the first 
(last) lattice site is put in contact with a fixed " 0" state 
on its right (left). 

1. The case a = 2.5 

The time evolution of p{t) has been analyzed for several 
values of v in the range [2.0, 3.0], starting from high den- 
sity equilibrium states. Close to the theoretical critical 
value, v c = 2.341486..., the dynamics has been sampled 
over a time lapse up to 600r, in order to obtain a re- 
liable identification of the critical point. Despite finite 
size effects are expected to affect significantly MC dy- 
namics in models with long-range interactions like the 
LPS model, we have been able to obtain the numerical 
estimate v c = 2.34 ±0.02, which agrees very well the the- 
oretical expectation. Moreover, for v > v c (i.e. T > T c ) 
we have also verified that p(t) initially evolves according 
to the dynamics p — ~j[e p — e] (see FigfTJ). 



Ap/Ax -i- 




_ 2 I i I i I i I i I i I i I 

-"0.06 -0.05 -0.04 -0.03 -0.02 -0.01 

[e-e(p)] 

FIG. 1: p = Ap/Ar versus e — e(p) for a = 2.5 and v — 
2.4, 2.45, 2.5, 2.53. The dashed line is the best-fit of the slope 
(—7 ~ 30), common to all values of v . It has been drawn to 
guide the eyes for appreciating the extension of the transient 
evolution, during which the four lines overlap. 

For larger values of time, p(t) decreases linearly in 
time, as predicted by Eq. (|6.86p . This confirms the pres- 
ence of the phenomenon of gelation, according to which 
the system reaches a denutared state in a finite time (see 
Fig©. 

We have not been able to check the quantitative agree- 
ment with the theoretical predictions reported in sub- 
section A. In particular, we have not found an effective 
criterion for locating the crossover between the initial 
and the asymptotic dynamics of p(t). It is expected to 
occur at the critical time t c (p(t c ) — p c ), but in MC sim- 
ulations the crossover region extends over a long time 
lapse, which typically begins before and extends far be- 
yond t c . We want to point out that such a scenario is 



o.s - 
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FIG. 2: The density of H-bonds p as a function of r for 
a = 2.5 and v = 2.4, 2.45, 2.5, 2.53 (from top to bottom). 
Note that a linear decrease in time is eventually approached, 
thus yielding denaturation after a finite time (gelation phe- 
nomenon) . 

not unexpected: finite-size and memory effects cannot 
allow for the identification of a single transition point 
between the two dynamical regimes. This analysis could 
be slightly improved by considering much larger lattices 
and longer simulations, while maintaining at least the 
same number of averages over initial conditions. We did 
not proceed along this line, because the limit of our com- 
putational capabilities has already been reached with the 
values adopted in the reported simulations. 

This scenario does not change also for v < v c (T < T c ) . 
We still find a qualitative agreement with the theoreti- 
cal predictions based on the order parameter dynamics. 
Specifically, for small times the dynamics is again ruled 
by p = -l[e P ~ e] (see Fig© . 




[e-e(p)] 



FIG. 3: p = Ap/Ar versus e — e(p) for a — 2.5 and v = 
2.2, 2.25, 2.3, 2.33. The dashed line is the best-fit of the slope 
(—7 py 30), common to all values of v . It has been drawn to 
guide the eyes for appreciating the extension of the transient 
evolution, during which the four lines overlap. 

Afterwards, one observes the decay to a finite asymp- 
totic value p € > p c (see FigHJ) . According to Eq. (|6.87|) 
such a decay is predicted to be exponential, with the de- 
cay rate r given in Eq. (|6.88p . 

Also in this case a clear quantitative verification is 
prevented by the long lapse of time characterizing the 
crossover between the transient and the asymptotic dy- 
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FIG. 4: The density of H-bonds p as a function of r for a = 
2.5 and v = 2.2, 2.25, 2.3, 2.33 (from top to bottom). Note 
that the decay to a finite value p e slows down significantly as 
v approaches v c — 2.34 ± 0.02. 



namics. 

Further peculiar behaviors of the microscopic dynam- 
ics emerge for v > v c and p(0) < p c . In this case any ini- 
tial condition cannot correspond to an equilibrium state. 
Accordingly, hysteretic effects determine a growth of pit) 
from p(0), until a value close to p c is approached. Then, 
p(t) starts to decrease until a linear-in-time decay is 
eventually approached, as predicted by Eq. (|6.85p (as an 
example, see FigfS]) . 




FIG. 5: The density of H-bonds p as a function of r for 
a = 2.5 and v — 2.6. The evolution has been obtained by 
averaging over randomly seeded initial conditions of initial 
density p(0) = 0.8 and 0.4 . In order to better appreciate the 
hysteretic effect associated with the first-order nature of the 
phase transition (lower curve), we have reported the evolution 
over a relatively short time scale. Note also that the initial 
conditions affect only the transient evolution: both curves 
converge to the same asymptotic dynamics, which eventually 
turns out to a linear-in-time decay, yielding the gelation phe- 
nomenon. 



well with the theoretical expectation v c = 3.6060508.... 
On the other hand, the crossover between the transient 
and the asymptotic dynamics is found to extend much 
more than in the case a = 2.5. For instance, even 
very close to v c the transient behavior ruled by the law 
p = — j[e p — e] lasts over a few units of r (see Fig|5J|. Such 
a linear region reduces the more v is far from v c 



Ap/AT 



[e -e(p)] 

FIG. 6: p = Ap/Ar versus e — e(p) for a — 1.5 and v = 
3.59, 3.60, 3.61. A linear dependence can be attributed to 
the first few units of r. 

However, the long-time dynamics allows to identify 
the occurrence of the gelation phenomenon for v > v c 
(T > T c ). In particular, p eventually exhibits a decay in 
time which is faster than linear. In order to exemplify 
this behavior in Fig[7]we show p{j) for some values of v 
much larger than v c . The asymptotic behavior is ruled by 
a decay faster than a power-law, but also slower than an 
exponential. This finding confirms the qualitative agree- 
ment with the predicitons of the order parameter dynam- 
ics. 




FIG. 7: The density of H-bonds p as a function of r for 
a = 1.5 and v = 3.8, 4.5, 6.0 (from top to bottom) . The 
double-logarithmic scale shows an aymptotic decay in time 
which is faster than linear. 



The phenomenon of gelation is recovered also in this 
case, as predicted on the basis of the order parameter 
equation. Nonetheless, a quantitative analysis is pre- 
vented for the above mentioned reasons. 
2. The case a = 1.5 

Numerical analysis predicts a critical value of the con- 
trol parameter v c = 3.60 ± 0.02. This result agrees quite 



We want to point out that we have reported the results 
for large values for v, because this choice allows to reduce 
the duration of the crossover from the transient regime. 
In fact, for the same values of v reported in Fig[6] the 
long-time behavior sets in for much larger values of r. 
Moreover, note that also in these cases finite size effects 
do not allow to reach a completely denaturated state (p = 
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0). In fact the long time dynamics exhibits fluctuations 
around a low-density regime which is not shown in FigJT] 

Also for v < v c (i.e. T < T c ) the transient dynamics 
p = — "f[e p — e] lasts for very short times and crosses over 
to a slow relaxation to the asymptotic non-zero equilib- 
rium value of the density, p e . As for the case a = 2.5, a 
quantitative study concerning the scaling analysis asso- 
ciated with the exponential decay rate cannot be worked 
out satisfactorily with our computational resources. 

Despite all the difficulties inherent MC simulations, we 
can conclude that it exhibits a remarkable qualitative 
agreement with the order parameter dynamics discussed 
in this Section. 



VII. SUMMARY 

In this paper we have derived new results for the clas- 
sical Lifson-Poland-Sheraga model of DNA denaturation. 
In the first part we have reviewed the main equilibrium 
properties and completed the list of known results by the 
computation of one- and two-point correlation functions 
and structure factors, including scaling formulas for the 
latter. In the second part we have investigated the dy- 
namical properties of this model, both numerically and 
analytically. The Langevin equation written for the den- 
sity of H-bonds has been solved with and without the 
noise term, in the second case in a linear approximation. 
Scaling laws for different values of a have been obtained, 
and a gelation phenomenon — arrival to equilibrium in 
a finite time — above the critical temperature has been 
identified. Remarkably, in the cases when comparison 
was possible, we have found that the predictions of the 
phenomenological theory were practically indistinguish- 
able from the numerical findings. 

Our results suggest that the LPS model belongs to a 



new universality class, in the language of critical phe- 
nomena, so that most of its details do not matter near 
the critical point. This is confirmed by both the static 
structure factor and the dynamics of the order parame- 
ter. Concerning the dynamics, the surprising success of 
the equation of the order parameter compared with the 
more detailed description by the master equation has its 
counterpart in usual treatments of critical phenomena. 
A theoretical explanation of this universality might be 
given by a renormalization group analysis of models of 
the helix-coil transition that are more microscopic than 
the LPS one. However, the best test of universality would 
be experimental. We hope that our work will stimulate 
such experiments, so that the very nature of the DNA de- 
naturation transition could finally be clarified. A better 
understanding of the gelation phenomenon characteriz- 
ing the high-temperature phase is also missing. Is this 
phenomenon basically related to the long-range nature of 
the effective interaction of the LPS model? To our knowl- 
edge no systematic investigation has been performed for 
identifying such a phenomenon in biomolcculcs. 

Note added: After finishing the analytic part of our 
work two papers (l7| , [l8l | have appeared on the preprint 
archive, dealing with the dynamics of the helix-coil tran- 
sition. Our approach is, however, sensibly different. 
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